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Abstract. We investigate the behaviour of randomly cross-linked (co)polymer blends using a combination 
of replica theory and large-scale molecular dynamics simulations. In particular, we derive the analogue of 
the random phase approximation for systems with quenched disorder and show how the required correla- 
tion functions can be calculated efficiently By post-processing simulation data for homopolymer networks 
we are able to describe neutron scattering measurements in heterogeneous systems without resorting to 
microscopic detail and otherwise unphysical assumptions. We obtain structure function data which illus- 
trate the expected microphase separation and contain system-specific information relating to the intrinsic 
length scales of our networks. 

PACS. 61.25.Hq polymer solutions - 61.43.j disordered solids - 64.75.+g phase separation - 05. 40. a 
fluctuation phenomena in statistical physics 
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I/*-) . 1 Introduction which heterogeneity does not play a role in the distribu- 

T^- 1 tion of cross-links. 

• ' Highly concentrated liquids comprising long polymeric chains 
)Q ! can undergo a process of cross-linking, resulting in the An f PP roach a ™ ed at circumventing the huge compu- 

formation of a disordered solid with a frozen memory tational expense of simulating blends subject to quenched 

1 of its preparation conditions. The cross-linking creates disorder involves devising a theoretical framework within 

' quenched connective and topological disorder which plays whl , cn heterogeneity may be incorporated after cross-linking, 

a role familiar from glassy systems in breaking the trans- wlthout rel y m S on the exhaustive computation required m 

lational symmetry of the initial liquid state. In the case brute-force simulations (see e.g. Ref 0). In this respect, 

of polymer blends, it prohibits macroscopic phase separa- a random P hase approximation (RPA) for the monomeric 

Cd tion so that demixing can only occur on a microscale . densit y Actuations is a convenient tool 0. However, one 

The resulting microstructure is intimately connected with nmst take , care to , modlf y the standard RPA expression 

the macroscopic elastic properties of the solid and proves often used to anal y se scattering experiments, before ap- 
essential to the correct interpretation of neutron scatter- :t to the s y stem wlth quenched disorder [H]. The 

ing experiments on dense cross-linked heterogeneous melts modifications arise from the fact that the density variables 

®. For example, partial deuteration (or 'labeling') is often fluctuate about some non-zero mean determined by the 

employed in scattering investigations of polymer blends [3, disorder, which effectively correlates fluctuations associ- 

0| . This process can result in weak but measurable effec- ated Wlth dln ? r , ent ™ VeCt ° rS ' nece ssitatmg an average 

tive interactions between labeled chain sections, inducing over q uenc hed disorder. 

a microstructure which influences the system's dynamics We construct an appropriate formalism, with the view 

and response to strain. to examine microphase separation due to interactions be- 

Simulations offer a powerful means by which to com- tween components in the melt, introduced after formation 

pare with the results of such experiments, claiming an ad- of the network. We utilise extensive foundations laid pre- 

ditional advantage in complete characterisation and con- viously for replica formalism in microscopic models [101 

trol of the microscopic state [5]. In particular, they are |TT| and field theories p~2l!13j . Whilst our theory follows 

capable of isolating (and thus negating) effects due to in- in the spirit of these studies, it does so without encod- 

homogeneities in the system prior to cross-linking. In this ing details about the microscopic structure of the chains, 

way, one can use simulations to investigate systems for This omission renders it applicable to our simulation data 
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which in turn provides insight into the effects of system- 
specific parameters on the outcome of (simulated) scat- 
tering experiments. Herein lies the strength of this work; 
we sample correlation functions for density fluctuations in 
homogeneous networks representing a broad range of net- 
work types and cross-linking densities [T41IT5] . From these 
we can infer the behaviour of heterogeneous cross-linked 
melts using our RPA formalism. In addition, we have de- 
veloped an inexpensive means of sampling the correlation 
functions in randomly labeled homogeneous systems. We 
can exploit this combination in order to address previously 
unresolved questions, such as the exact behaviour of the 
structure factor near microphase separation [T^flo]. 

In the next section, we define the system simulated via 
the molecular dynamics of cross-linked polymer systems 
and describe in detail the method by which we take care 
of the quenched disorder. This section is followed by an 
application of spin-glass replica formalism [17] to our data 
in order to describe neutron scattering measurements in 
heterogeneous systems. We conclude with a discussion of 
microphase behaviour and its implications. 



2 Simulated scattering from homogeneous 
networks 

We study data from extensive molecular dynamics simu- 
lations of randomly cross-linked bead-spring polymer net- 
works [13]. For completeness, we briefly characterise the 
system in an Appendix. Partial deuteration of the net- 
works is implemented by randomly distributing labels of 
variable length, Ng, on precursor chains. Of the total num- 
ber of all available labels, N a , we label a fraction, 0, such 
that the system contains <j)N a labels in all. In this way, we 
can assume there are no correlations between labeled and 
unlabeled segments. This plays an essential role in sim- 
plifying our analysis and renders the calculation of corre- 
lation functions far more straightforward than the anal- 
ogous process for an end- linked triblock copolymer melt, 
for example. 

We obtain correlation functions corresponding to mea- 
surements from scattering experiments in terms of label 
density, p(r), which is defined by the position of monomer 
I on label a such that, in Fourier space, 

a £ 

The scattering observed in the precursor melt is then de- 
scribed by the structure function, 

5 (q,k) = -i-(p q ^) , (2) 

where (. . . } represents a thermal average over annealed 
variables in the context of the molecular dynamics simu- 
lations and the subscript '0' reminds us that the system 
is homogeneous. 



The quenched disorder creates correlations between 
density fluctuations associated with different wave vec- 
tors. The translational invariance of the simulated precur- 
sor state is lost upon cross-linking and the density vari- 
ables fluctuate about some non-zero mean determined by 
the disorder, 

( Pq ) ^0 for q^O, 

(pqPk)o ¥= ° for q^ k - 

An average over disorder restores the invariance expected 
in macroscopic systems, and resulting correlation func- 
tions may be expressed as, 



<Pq) = °. ( 3a ) 

(PqPk) ° = N a mS {q)5^, (3b) 

(Pq)o (Pt)o° = N a Nir (q)S^ (3c) 

for q ^ 0, where the overline denotes a quenched disorder 
average. These equations describe respectively the vanish- 
ing mean of density fluctuations, the structure factor for 
the homogeneous system and the overlap function familiar 
from spin glass physics [18] . The latter term in particu- 
lar arises from the broken translational symmetry in the 
quenched system. 

A straightforward way to calculate these averages would 
be to (randomly) assign fixed labels to chains and to aver- 
age over a (suitably large) number of independent realisa- 
tions of the vulcanisation process. In this case, we would 
need sum over the fraction of labeled monomers for all 
realisations, 

y \aa' \ W / / y 

(4a) 

Y-l X MN a I Nt \ I N e \ \ 

y \aa' \ I /q \ V I Q/ ., 

(4b) ' 

where y indexes disorder realisations and the correlation 
functions are calculated within each realisation. One way 
to circumvent excessive summing involves supposing that 
we can obtain equivalent results by taking a single distri- 
bution of cross-links and summing correlation functions 
over different ways of labeling the system, 

X /<f>N a \ ("nJ <pN a 

£<->o ^jh E E <->o- (5) 

y \ao' Jy \(j>N a ) x act' Ex 

where x enumerates ways of choosing labels. The indices a 
and a' are then taken from the set specified by x. In terms 
of scattering measurements we should be able to sum over 
exactly the same functions this way, given a large enough 
system. Nevertheless, we are still obliged to sum up over 
the fraction of labeled monomers as many times as there 
are ways of choosing those labels and this is unnecessar- 
ily computationally exhaustive. Fortunately this method 
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leads to a more efficient means of calculating these aver- 
ages. The sum over different ways to choose 4>N a labels can 
be rewritten as a weighted sum over same- and different- 
chain correlations between all possible N a labels, 



average, we would be wrongly assuming that the quenched 
variables (namely, the cross-links) fluctuate on the same 
timescale as the density variables. Instead we rewrite Eq. (JT]) 



S (q)<5 q k = 



N a N t 



a \ W 



a a'^a \ W 



rf— k r 




(6a) 



-To(q)<5 q k = 



N a N t 

N a N a i I N, 



E(E eiqrf / (E e ~ 4kr 

a \ £ /o \ £' 



4>N a 1 \ - \ - / 

"ivw E E (E* 



iqr f 



E' 



-ik-r™, 




where the primed sum indicates that each t and £' do not 
appear on the same chain. This drastically simplifies com- 
putation, since one need only consider a single realisation 
of connective and topological disorder. The calculation of 
all possible pair correlations within this realisation yields 
the disorder-averaged correlations appearing in Eqs. ([3]). 



3 Replica Random Phase Approximation 

We now construct a formalism appropriate for calculating 
disorder-averaged observables of heterogeneous systems, 
using those obtained from simulating homogeneous sys- 
tems. As in other systems characterised by quenched dis- 
order, we face the task of averaging over the logarithm 
of a disorder-dependent partition function. To accomplish 
this, we employ the techniques and language of replicas 
[18] , well-known from spin glass theory [19] . This approach 
has been used extensively in microscopic formulations of 
analogous polymer systems and is particularly convenient 
in our case [TOlfTT] . 

The problem amounts to one of taking the disorder 
average of a given observable, (^4) - This quantity is itself 
thermally averaged via an integral over fluctuating density 
variables, with Boltzmann weight determined by a Hamil- 
tonian, Ho[/Oq], such that the disorder average takes the 
form, 



where, 



{A} =V / P Pq A[p q ]exp(-/W [ Pq ]) 



Z = / P Pq exp(-/W [p q ])- 



(7) 



(8) 



The implicit dependence of the normalising partition func- 
tion, Z , on the distribution of quenched cross- links poses 
an analytical challenge. By removing it from the disorder 



(A) = lim ZJ- 1 / P Pq A[p q ]exp(-/3WoM 

n->0 J 

and make the following identification, 

z o = J (fl v P a ^j ex P (-0 E n ° K^j ■ 



(9) 



(10) 



This is the essence of the replica trick and may be inter- 
preted as an n-fold replication of the system. It should 
be stressed that every Hamiltonian appearing in Eq. (|10p 
corresponds to the same realisation of disorder. We can 
now write, 



!™ / (j[ vp ^j ^Ki^pf^E^K]) . 

(11) 

where the choice of replica label in the functional argu- 
ment of A is arbitrary. We further simplify matters by 
introducing an effective Hamiltonian given by, 



Hf[pl] = -f3- l \n e -"S:WoI^l 



(12) 



which describes the energy of the disorder-averaged sys- 
tem. Now the average over observable A can be recast in 
the form, 



(A 




Il^q A K1 exp(-/3^ ff [p q ]). (13) 



Note that the Hamiltonian Ho [p q ] for a given realisation of 
quenched disorder is replica diagonal, whilst the effective 
Hamiltonian Ho*!/ 9 ^] retains interactions between differ- 
ent replicas. 

We can use this recipe for weighted averages to de- 
scribe the density fluctuations calculated in our simula- 
tions after a disorder average. In this parametrisation, 
Eqs. ([3]) can be rewritten as, 



\eff 



{p a ^ a T ^N a N e S (cD, 

>cff 



/ a *b \ u 



N a N e r (q), a + b. 



(14a) 
(14b) 
(14c) 



for q 7^ 0. We need only choose a reasonable Ansatz for 
7i§ ff [p q ] in order to recover these correlation functions us- 
ing Eq. ([T3]) . Invoking the RPA, we assume Gaussian fluc- 
tuations in density. Under a further assumption of symme- 
try under permutation of replicas, we can uncouple wave 
vector space by coupling replicas to write, 



1 



wf[p q ] = iEE' o q K «f( t iK' 



(15) 



ah 
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with [K- X (q)]ab - (p^Tq such that 
[K-^q)]^ = N a N e (5 ab S (q) + (1 - 



The inverse corresponding to Eq. 
for the kernel in the form, 



~S ab )r (q)). (16) 
161) yields an expression 



K afc (q) 
(1- 



5 (q) + (n-2)r (q) 



N a N e (5 (q) 

Sab) 



r (q))(5 (q) + (n 



l)A(q)) 



JV Q JV* (5 (q) - r (q))(5 (q) + (n - l)r (q)) ' 

(17) 

In order to investigate a heterogeneous analogue of 



The structure of the equations may be directly com- 
pared with results from microscopic descriptions of end- 
linked [16] and cross- linked [13] systems. The advantage 
with this formulation comes from the fact that we are 
free to explore material parameter space and deduce the 
behaviour of the system in different limits provided by 
the simulations. The derived equations clearly show the 
change effected in the behaviour of the system upon the 
introduction of inhomogeneities. The two heterogeneous 
correlation functions are formed by shifting their homo- 
geneous counterparts by a self-energy term controlled by 
the interaction potential 
form in Fourier space. 



Vq, which retains its functional 



this system, we now introduce an interaction potential, Vq 4 Results 
within the formalism. This inclusion is straightforward in 
the case for which the potential is a function of monomer 
separation only (and therefore diagonal in wave vector 
space). We perform the same formal treatment as for the 
homogeneous system, by first introducing the replica trick 
and then averaging with respect to disorder. The limit of 
vanishing replicas is trivial and physically reasonable at 
this stage of the calculation and in taking it we obtain the 
effective Hamiltonian, 



Hf[, 



1 



=^m + ^EEw, as) 



2Q 

q a— 1 

where the system is enclosed within volume £2 7 and 
-Vtq)]-, 1 =N a N e (6abSv(<l) + (l 



K 



(19) 

with [V(q)] a f, = S a bVq/Q. Using the result from the ho- 
mogeneous system in Eq. (|17[) , we can rewrite Eq. (|19p in 
the form, 



[K(q) 
5o(q) 



•VCq)]^ 
h gVq(S (q) 



l + gVq(S (q) 



Io(q)))(l + gVq(So(q) + (n 
N a N e (l - 6 ab )r (q) 



1 + gVq(S (q) - r (q)))(l 



where g = N a N^/f2 defines the density of the system, and 
from this expression we can infer that, 



Sv(q) = 



5 (q)-r (q) 



1 + gVq(S Q (q) 



(1 



1 



gVq(S (q) 

in the replica limit where n — > 0. 



The simulations provide data for several different param- 
eter values and we can further tune our output using pa- 
rameters arising from the post-processing. We begin by 
considering systems with label length Nt = 100 and vary- 
ing the number of cross-links, which in turn determines 
their average chemical separation, N x . This variable in 
particular plays an important role in the behaviour of the 
overlap function, -To(q). In Fig. QJa), we show data for 
the overlap function for different values of N x . Data for 
iSo(q) is inset here for comparison. We are restricted to a 
finite simulation box and thus cannot capture behaviour 
in the very long wavelength limit. However, we note the 
convergence of both correlation functions in this limit to 
finite values which we normalise to unity. It should be em- 
phasised that different sets of <5>o(q) data superimpose one 
another, implying that the data show no A^-dependence. 
It is clear from the r (q) data that the overlap is depen- 
dent on N x , at least below a threshold of the order of the 
melt entanglement length of L e = 70 20J . 

We use the expression given in Eq. (|21a[) to post-process 
the data and obtain structure function data for the het- 
erogeneous system. For local attractive interactions be- 
tween labeled monomers, we set Vq = — x and observe a 
divergence in <Sy(q) at finite wave vector where the dif- 
ference 5o(q) — ^o(q) forms a maximum. These maxima 
are illustrated in Fig. [ljb), for different values of N x . In 
Fig. [2] we show the post-processed data for different pa- 
, rameter values. Once again we note the convergence to a 
gVq(So(q) + (n — l)Jo(q)) n on- vanishing scattering intensity in line with early exper- 

( 20 mients [2] and at odds with de Gennes original prediction 
PP. This discrepancy has been the subject of much inves- 
tigation [21,22,16,6,12,13 , but the long wavelength ten- 
dency in our simulations is unambiguously non- vanishing. 

In all plots in Fig. [2] the emergence of a character- 
istic length scale represented by a peak at finite q may 
be interpreted as the establishment of pronounced spatial 
correlations between labeled chain densities. The position 
and size of peaks are strongly dependent on our position in 
parameter space. In each plot we vary a single parameter 
(interaction strength x, label fraction <fi and average cross- 
link separation N x respectively) and examine the wave 
vector dependence of data leading towards divergence at 
characteristic values qs and \s = (<->o(qs) — A(qs)) -1 - 



N a N e 8 ab x 
- r Q (q))(S (q) + (n 



l)A(q)) 



l)A(q)) 



^o(q)) 

r (q 



■ gVq{S (q) 

r (q) 



r (q)W 



r (qW 



(21a) 



(21b) 



A. V. Klopper et al.: Microphase separation in cross-linked polymer blends 



5 




Fig. 1. (a) Normalised structure function <5>o(q) (inset) and overlap function -To(q) data from molecular dynamics simulations 
of systems with different average cross-link separation iV^. Here, a is the length scale defined by the Lennard- Jones potential 
used in the simulations (see Appendix for details). Overlap data show strong dependence on N x when it is shorter than the 
entanglement length, L e = 70, whereas those for weakly cross- linked systems, as well as all structure function data, show no such 
dependence. Solid lines in main plot interpolate data. Solid line in inset plot corresponds to the Debye function, (b) Difference 
between normalised correlation functions >So(q) — ^o(q) peaks at characteristic wave vector qs, which is dependent on cross- link 
separation, N x . 



In each case we also plot the value of qs, which is ob- 
tained by fitting the appropriate Sq (q) — Jb (q) curve and 
locating a maximum. 

In Fig. [Ula) we see the effect of increasing \ f° r a given 
fraction of labels and cross-links. There is a pronunciation 
of the peak leading to divergence and a corresponding shift 
in the wave vector value at which this peak occurs. A sim- 
ilar trend can be seen as the fraction of labels increases, 
as shown in Fig. [2h). In this case, we see a characteris- 
tic length scale evident in systems with 20% of polymers 
labeled for a given interaction strength. Under the same 
conditions, the structure function for the corresponding 
system with 5% labeling shows no such tendency and in- 
stead the data resemble those of the homogeneous system. 
This disparity is more pronounced when comparing sys- 
tems with different average cross-link separation, as seen 
in Fig. [21c). We do not detect a divergence in strongly 
cross-linked systems, even as identical weakly cross-linked 
systems show a marked effect. Note that the inversion of 
data for N x = 100 and N x = 200 in Fig. Hfc) should be 
interpreted as an artifact of the simulation statistics. 



5 Discussion and conclusions 

We note that qs is independent of \ an d 4>i the latter 
implying that the length scale describing a single label 
interacting with itself is indistinguishable from that de- 
scribing several labels interacting. This in turn suggests 
that the divergence of the structure function is in fact a 
single-label effect. This observation is supported by other 
models but experimental and simulation results give the 
impression that more cooperative effects are prevalent in 
such systems [5], However, we also note that the approach 
to separation is (^-dependent, as shown in Fig. Hfb). This 
may be interpreted in the following way. The system is 



prepared in the initial homogeneous melt state and cross- 
links arc formed. Simulated scattering measurements are 
indistinguishable from those of a polymeric liquid. Inter- 
actions are introduced between labels and the structure 
function is deformed slightly, but the physical properties 
remain largely unchanged. As interactions increase, some 
structure emerges at an intermediate length scale corre- 
sponding to a finite wave vector. Interactions are increased 
still further and labeled monomers become more attracted 
to one another, contracting the characteristic length scale 
and thus shifting the value of the corresponding wave vec- 
tor. When the interaction strength reaches its spinodal 
value, the length scale of 'domains' coincides with the in- 
trinsic length scale of the system, related to the length 
scale of single-label collapse. The value of qs is indepen- 
dent of labeling fraction because the self- interacting length 
scale dominates the shape of the structure function. 

Conversely, the value of qs is strongly dependent on 
N x when N x < L e = 70 as illustrated by Figs. [He) and[3] 
In the latter, the dependence of qs on N x appears to obey 
an inverse square root law as seen in experiment [2] and as 
predicted by classical rubber elasticity theory [1] . This be- 
haviour holds for different label lengths and breaks down 
for systems in which the average cross-link separation ex- 
ceeds the entanglement length of the system, as indicated 
on the plot. We note in Fig. [4] that these data collapse 
when scaled with label length-dependent radii of gyration, 
(Rg) {Nji), and average cross-link confinement, dx [141115) . 

Finally, we find that we can use our post-processing 
formalism to determine exactly the character of the di- 
vergence for x near \S- The mean-field result predict- 
ing ixs — x) behaviour was shown to describe exper- 
imental data well [2], but a (xs — x)~ 2 divergence has 
been formally derived [TBIUS] , We are able to locate the 
crossover from mean-field to higher order behaviour in the 
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Fig. 2. Post-processed structure functions for systems with (a) 
= 0.1, N x = 100, variable X ; (b) X = 0.4, iV* = 100, variable 
0; and (c) = 0.1, x = 0.7, variable N x . Broken lines indicate 
the value of qs for respective parameters. 



parameter space of our simulations. In analogy with pre- 
vious studies [T|l2|ll6j. we rewrite 5y(qs) in terms of the 
strength of the attractive interaction, %■> as j 



5y(qs,x) = 



1 



| x 2 s r Q ( qs ) 

Xs - X ^ (xs - X) 2 



(22) 



0.75 



5S< 

b 




0.25 



Fig. 3. Dependence of spinodal wave vector qs on average 
cross-link separation N x . Solid line fits Nt = 100 data; bro- 
ken line indicates length of entanglement, L e . Again, a is the 
length scale defined by the Lennard- Jones potential used in 
simulations (see Appendix for details). 



2.5 ~ 

2 - 

1.5 - 



0.5 - 



□ N x = 10 

O N< = 2Q 

A N x = 100 

jV, = 200 



10 

d x 2 {Rb(N ( ) 



100 



Fig. 4. Data collapse when scaled with A^-dependent radii of 
gyration, (-Rg) (Nt), and average cross-link confinement, dx- 



We define a parameter measuring 'distance' from the di- 
vergent interaction strength, e = 1 1 — x/Xs|j an d note that 
the crossover occurs at some e* given by, 



1 - 



r 

Xs 



xs-fb(qs)- 



(23) 



However, data from our simulations give e* S [0.66, 0.74] 
for different values of N x , indicating that higher order be- 
haviour dominates our systems near the divergence. Ex- 
perimental data [2] showing a linear relationship between 
5y(qs) _1 and an interaction parameter analogous to x 
may not depict the truly divergent region of the systems' 
phase space, as illustrated by Fig.[5l and noted in Ref. [16]. 
Here, we show that the linear relationship breaks down as 
X approaches xs- One may compare these curves with 
those in Fig. 7 of Ref. [2] . Our control over network char- 
acteristics via simulation proves particularly advantageous 
here; one can clearly see that the crossover to higher or- 
der behaviour occurs over a broader range of x values in 
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Fig. 5. Crossover from (\s — x) _1 to (xs — x)~ 2 behaviour 
occurs as x approaches xs, indicated here by broken lines. A 
linear curve is plotted for comparison. 



strongly cross-linked systems. This change is much harder 
to predict in more melt-like cases for which the average 
strand length exceeds that of entanglement. 

To conclude, we have developed a means of examining 
microphase separation in polymer systems with quenched 
disorder without the need for exhaustive simulation. The 
formalism we use introduces heterogeneity in a controlled 
way and emerging results show a reproduction of exper- 
imental data. The observation of microphase separation 
in systems interacting via a localised potential calls for 
further analysis using these computational and analytical 
methods, particularly with respect to the accessible pa- 
rameter space of the simulations. The preservation of the 
functional form of the interaction potential suggests the 
formalism may even be applicable to ferrogels [23j and 
related systems involving non-local interactions. 



Table 1. Simulated network characterisation in terms of tar- 
get strand length, N x ; effective strand length, (N 3 ); effective 
functionality of cross-links, (/); defect-free strand percentage, 
d.f.s.; and defect-free bead percentage, d.f.b., where dangling 
ends constitute a network defect. 



N x 


(Ns) 


(/> 


d.f.s. 


d.f.b. 


10 


13 


5.1 


99% 


99% 


20 


23 


4.4 


99% 


99% 


100 


101 


4.1 


94% 


93% 


200 


192 


4.0 


89% 


87% 



for r < 2 1/6 cr. 



The physical units for distance, energy and time are 
defined by parameters a, e and r = o\J m/e respectively 
and all beads have mass m = 1. The standard Kremer- 
Grest force field [2~4"1T2"5] parameters have values Rq — 1.5a 
and k = 30eer~ 2 . These parameter choices yield a high en- 
ergy barrier for chain crossing which ensures conservation 
of the topological melt state. The Langevin dynamics is 
integrated using the LAMMPS code [26] with tempera- 
ture ksT — le, monomeric friction r = 0.5r _1 and time 
step At = 0.012t. The resulting model polymers map onto 
those found in natural rubber [27] . 

Melts of 80 chains each with 3500 beads are generated 
with bead density g = 0.85a-" 3 and equilibrated as de- 
scribed in Ref. [28]. The melts are instantaneously cross- 
linked by adding between 700 and 14000 FENE bonds 
between random pairs of beads separated by less than 
the reaction radius of 1.3<r. This corresponds to nominal 
strand lengths N x between 10 and 200 monomers. The re- 
sulting randomly cross-linked networks are characterised 
by an exponential strand length distribution [14] . Table Q] 
summarises our characterisation of the network. 
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A Details of the simulations 

Polymer chains are modeled as sequences of beads con- 
nected by FENE springs, while excluded volume inter- 
actions between beads are described by the shifted and 
truncated 12-6 Lennard- Jones potential, 

c/ F EN E (r) = — Y ln y~Ri) for r<i?0 ' (24) 
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